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ABSTRACT 

Input-output, growth-decay, production-consumption type situations abound in many practical prob- 
lems. When the input and output variables are independently gamma distributed, various aspects of the 
residual effect are already tackled by the author. Matrix-variate analogues, their connections to quadratic 
and bilinear forms, matrix-variate Whittaker functions, and many properties of such matrix-variate functions, 
in the real as well as complex case, are discussed here. 



1. INTRODUCTION 



If the input variable x\ and the output variable xi are independently gamma distributed then it is 
known that the residual effect y = x\ — X2 has a density which can be expressed in terms of Whittaker 
functions in the scalar case. A particular case is a Laplace density. When Xj has the density 



fj( x i) 



xj > 0, Re(aj) > 0, Re{/3j) > 0, j = 1,2 



(1.1) 



and fj(xj) = elsewhere, where Re(-) denotes the real part of (•), the density of y = x\ — X2, denoted by 
f{y), is given by 
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and W, .(•) is a Whittaker function, see also Mathai (1993a). 
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Independent sums of such residual effects are shown to be connected to the distributions of quadratic 
and bilinear forms in Gaussian random variables, see Mathai and Provost (1992), Mathai, Provost and 
Hayakawa (1995) and Rao (1973). Laplacianness of quadratic and bilinear forms or the conditions under 
which a quadratic or bilinear form is distributed as a Laplace variable or as a gamma difference is also 
investigated. Let Sbeapxp non-singular central Wishart matrix and let it be partitioned as 



S = 



Su S12 

S21 S22 



(1.3) 
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where Sn is r x r, r < p then it is shown that the distribution of Si 2 is associated with that of certain 
bilinear forms. Note that S\ 2 can also be looked upon as a certain covariance structure when S is a corrected 
sample sum of products matrix. A large number of results in this category may be seen from Mathai (1993b) 
and Mathai, Provost and Hayakawa (1995). 

Solar neutrinos are captured on earth and the data for over twenty years are available. The summarized 
data, by taking a five point moving average to smooth out local disturbances, show a peculiar pattern of a 
9-year cyclic behavior where within each cycle the pattern is a slow rising and rapidly decreasing curve with 
one real peak and several small humps in between. This type of pattern is also seen in many other problems 
such as the production of melatonin in human body. In Haubold and Mathai (1995) it is shown that this 
type of curve could be generated by combining three residual variables where the input and output variables 
are of gamma type, and some heuristic interpretations of these variables are also given there. 

One matrix-variate analogue of gamma type input and gamma type output gives the residual variable 
Y = Xi — X 2 where Xi and X 2 are independently distributed p x p matrix-variate real gamma random 
variables with the densities 

#j) = J^lA-r-^e-W^, (1.4) 
for Xj = Xj > 0, Bj = Bj > 0, Re(aj) > and gj(Xj) = 0, j = 1, 2 elsewhere, where for example, 

T» = ^-^ 4 T(a)T (a ...T (a , Re(a) > (1.5) 

is the matrix-variate gamma function in the real case, with the integral representation 

T p {a) = f \X\ a -^e- tr{x) AX. (1.6) 

JX=X'>0 

When Bj = ^Aj for some Aj = Aj > one has the non-singular central Wishart density from (1.4). Then 
if Xi and X 2 are two independent central Wishart matrices then the residual variable is 

Y = X 1 -X 2 . 

One can look at element- wise differences in X\ — X 2 and consider the configuration of residual variables or 
one can consider the difference in terms of the definiteness of the p x p matrices X\ and X 2 . Since X\ and 
X 2 are real symmetric positive definite matrices in (1.4), Y could be positive definite or negative definite or 
indefinite or scmidefinite. A special situation will be when X\ and X 2 are oriented matrices such that Y 
is either nonnegative definite or —Y is nonncgative definite. In this case one can proceed in the following 
way to compute the density of Y . The joint density of X\ and X 2 , denoted by f(X\, X 2 ), is the product of 
fi(Xi) and f 2 (X 2 ) due to statistical independence. That is, 

f(X 1 ,X 2 )dX 1 dX 2 = S\X 1 \ ai - E i 1 \X 2 \ a2 ' E i 1 e^ BlXl+B2X2 Mx 1 dX 2 (1.7) 

where 

{ = jW.. (L8) 

rp(ai)r p (a 2 ) 

Consider the transformation U = X\ — X 2 and V — X\. The Jacobian is unity and the joint density of U 
and V, denoted by g(U, V), is given by 

g(U, V)dUdV = SIV] 011 -^ \V - Ul^-^e-^^+^-^dUdV. (1.9) 

The density of U, denoted by g(U), is available by integrating out V from (1.9). Note that when U is oriented 
as described earlier there are only two possible situations, U > 0, V > U and U < 0, V > 0. Let us consider 
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these two situations separately. For U > 0, V > U let 



gi(U)= f g(U,V)dV 
Jv>u 

= S [ |y|«i-4 i |y_ C /|«2-4 i e- tr (s 1 ^+B 2 (v-c/)) d y 
Jv>u 

= S p\ ai+a2 -2±± e -tr(B 1 U) f \ Z \*2-^ 

Jz>a 

„+i ~tr(ui(B 1 +B 2 )uiz) 

I + Z ai ~~e V A 



x |/ + Z| Q1 - — c V " 7dZ (1.10) 

by making the transformations Y = V — U and Z = U~?YU~^ for fixed U. The integral in (1.10) even when 
p = 1, that is in the scalar case, is difficult. In the scalar case it can be evaluated in terms of a Whittaker 
function which can then be represented as a linear combination of two confluent hypcrgeometric functions 
at least in some special cases. But when p > 1 no such representation is available due to the fact that when 
p > 1 

/ "/ + / 

JX=X'>0 J0<X=X'<1 Jx=X'>I 

because < X < I implies only that the eigenvalues of X are between zero and one and X > I gives the 
eigenvalues greater than or equal to one. But a whole class of matrices with some eigenvalues less than one 
and others greater than or equal to one are left out from the set X = X' > 0. Hence the integral in (1.10) 
is not available even from a matrix-variate version of a confluent hypergeometric function. We need the 
Whittaker function of matrix argument to deal with this integral. This will be defined next. 



2. WHITTAKER FUNCTION OF MATRIX ARGUMENT 



When studying some properties of the matrix-variate Laplace transform, as applied to special functions 
of matrix argument, it was necessary to come up with an expression for an integral of the following type. 

/ |AT-4V tr < x )dA = f \Y + A\ a -'?e- tt WdY 

JX>A JY>0 

= \A\ a c-^ A ) f \I + Z\ a - P -^e^ AZ Mz 
Jz>o 

by making the transformations Y = X — A, Z = A~^YA~?. Hence Mathai and Pederzoli (1996) defined a 
function M(a, (3; A) as follows: 

M(a,0;A)= f lAT^ \I + Xf^ e~^ Ax ^ dX (2.1) 

Jx=X'>0 

for A = A' > 0, Re(a) > and no restrictions on {3. Then we may observe the following from the definition 
itself. 

M (a,?±±;A) = \A\- a T p (a), Re(a) > (2.2) 

But for p = 1, that is , in the scalar variable case the function M {■,■;■) is associated with a Whittaker 
function. Hence we define a Whittaker function in terms of M(-, •; •) as follows: 

M(n,v;A) = lAI-^r.^e^)^^^^.^^^) (2.3) 
for A = A' > 0, Rc(^i) > and no restriction on v. In terms of the integral representation we have 

/ \Z\'"- E ^-\I + Z\ v - 2 ^-e- t ^ Az ^dZ 
Jz>o 

= |Ar^r p ( Ai )e^ r ^IU, ( _ Al)! , (iy+M _ (p+1)/2) (yl). (2.4) 
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By using (2.4) one can establish a number of results on Whittaker function which will generalize the cor- 
responding univariate results to the matrix-variate case. For the sake of illustration one such result will be 
listed here, before going back to the density in (1.10). This result follows from the above definition itself. 



Theorem 2.1. For Re(/3 - a) > = A' > 0, 

f |z|/3 _ a _41 |/+z|a+ 0_£+l e _ tr(A z) dz 

Jz>o 

= \A\-f~^T p (/? - a + £±i) eW% aJjW 

This follows from the definition itself by taking a = \{y — /j,) and (3 = | (y + /j, — ^-)- 
Compare the integral in (1.10) with the result in Theorem 2.1 to obtain 

9l (U) = S\B, + B 2 \-^T p {a 2 )h 1 {U) 
where 

h,{U) =e^ r ^- Bl ) c/ )|C/|^- £ * i W |(Q1 _ a2)ii(ai+a2 _ (p+1)/2) ((B 1+ B 2 )^U(B 1 + B 2 )^ .(2.5) 
Now for U < 0, V > let 

g 2 (U)= f g(U,V)dV. 
Jv>o 

Write W = -U so that W = W > 0. After taking out make the change Z = W~iVW~i for fixed W 
to obtain 

g 2 {U) = s\w\ ai+a '-'^e-' t '^ vr > 

f p+ i p+1 -tr( wi(B 1 + B 2 )wiz ) 

Now in the light of Theorem 2.1, we have, for W = —U. 
92{U) = 6\Bi + B 2 \-^T p ( ai )h 2 (W) 

where 

h 2 (W) = iTTl^-^ei^^-^)^^^^^^^^^/^ ((Si + S 2 )^(B! + B 2 )*) . (2.6) 
Hence the density of Y = X\ — X 2 , in the oriented case discussed in (1.10), is given by the following 
Theorem 2.2. 

( c-ifniY), Y>0 
g (Y) = I c-ifi^W), W = -Y > (2.7) 
[ 0, elsewhere, 

where h\(Y) and h 2 (W) are given in (2.5) and (2.6) respectively, c = c\ + c 2 , with 

ci= / /n(y)dy and c 2 = / ^(wjdw. 

The density in (2.7) generalizes the univariate case Mathai-1,4 of Mathai (1993a, p. 34). 

Integrals of the type in (1.10) appear in a large variety of statistical distribution problems connectd 
with sum, difference and linear functions of positve random variables, especially of the gamma type, or when 
dealing with the upper part of the incomplete gamma function. The corresponding matrix-variate analogues 
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can be handled by using the definition of Whittaker function given above. When dealing with the distribution 
of a linear function of three independent real gamma variables one has to evaluate a double integral where 
the kernel part is of the type in (1.10). This leads to integrals involving Whittaker functions, see Mathai 
and Pedcrzoli (1996). (Distributions of linear functions of gamma variables are discussed in Mathai and 
Provost (1992)). Several results on integrals involving Whittaker functions of matrix argument are given by 
the author and his co-workers recently. For the sake of illustration one such result will be given here as a 
theorem. 



Theorem 2.3. For A = A' > 0, B = B' > 0, Re (7 ± (3) > 2^, Re ((3 - a) > 2f5, Z = Z' > 0, 



f \Z[<-^ L G- t " {A ^W a ,i3{BZ)AZ = 

JZ>Q 

r„ (7 + /? + ^)r p (7-/3 + ^) 
r p ( 7 - a + £±i) 



A+-B 



-(7+/3+4 1 ) 



1 



B^ 



for0<B2[A+-B) B2<Ior2B~2AB-2>I 



(2.8) 



where 



. R |-( 7+g+ *ji) r P + 7 + /3) r P (^ + 7 - /?) 
1 ' r P (H±i + 7- a ) 

x a F 1 ('^+7-/3,^+7 + /3;^+7-a;-C I 

for \\C\\ < 1 



C = B-2AB-z - -I 
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and ||(-)|| denotes a norm of (•). ,4 sufficient condition is that 



< B~^AB~i - < / or < - B~^AB~i < I. 



(2.9) 



Proof. Replace W ai p(BX) by its integral representation to obtain 

left side = I \Xp- E i 1 e~ tl{AX) W a ^{BX)dX 



|X|7+/3+4 i - £ J i 



-tr(AX)-ltr(BX) 



L>o r p (2±i+/J-a)' 

/ |Z|' 3 -«-4 i I J + Z |a+/3-H^ e -tr(BXZ) dzdA -_ 

iz>o 



The X-integral can be evaluated by using a gamma integral to get 



/ 

Jx>o 



|x|7+/3+ £+i_£+i e _ tr 



[((A+lB) + BZ)x} dx 



r P 7 + /J- 



p + 1 



for Rc( 7 + /3) > 



p-3 



BZ 



-(7+/3+4 1 ) 



(2.10) 



The ^-integral, denoted by /, is given by 



/ = / \zf- a -^\i + z\ 

Jz>o 



A + -B ) + BZ 



-(7+/3+4 1 ) 



dZ 



A+-B 



-( 7 +/5+£+l) 



/ \z\P- a - p -^\i + zr : V 

Jz>o 



I + B* [ A+-B) B*Z 



-(7+/3+4 1 ) 



dZ. 



Note that in the determinant we can also replace (A + \B) B by B^ {A + \B) si Write Z = (I - 

U)-iu(I-U)-i => I+Z- 1 = U- 1 |Z|-(P +1 )dZ = |J7|-(f +1 )dC/ ^ dZ = |7-[/|-(P +1 )dt/ and < U < I. 
Then 



/ 



A+-B 



-(7+/3+4 1 



\u\ 



0-OL- 



0<U<I 



I-B* [A+-B] B^\U 



■| 7 _ j7|7-/»- 
-(7+/3+4 1 ) 



dC/. 



Now we can write the integral as a 2-F1 by using Mathai (1993a, p. 179). Then 

f = \A + ±B\-(n + P +t i L ) 

Tp ((3- a +^)r p ( 7 -/3+2±l) 



T p (7 - a H 
2 F, ( P±l+p- a , 



-H±l) 
p+1 



+ 7 + /?; 



for < 7- Bi (A+ \B) 1 fii < 7 => < (A + ±B) 1 < 7 2B-3 AB~3 > 7. Substituting back, 
one gamma is cancelled, we get (2.8). But observe that 



A+-B ) +BZ 



A--B)+B(I-Z) 



= \ B \ I 1 - z \ 



I + - ^ (7 + zr 1 



Now 



/= | B |-(7+/3+^) /" IZI^-^IJ + Z^ 
iz>o 



7+ (s^AB - * - -J J (J + Z) 1 



^(7+/3+4 i ) 



dZ. 



Put £7 = (7 + Z)- 1 =>dZ = |C/|-(P +1 )d?7, \Z\ = \U\~ l \I - U\, |7 + Z\ = \U\-\0 <U < I. Then 



/ = | jB |-(7+/3+4 i ) /" 

Jo 



7-/3- 



0<(7</ 



17 - |j + cu\~ (^+4 i ) 



where 

C = B-iAB-i - -I. 

2 

Evaluating this integral by using Mathai (1993a, p. 179) one has 

, _ , R |-( 7+ g + *P) r p + 7 - 0) T P + P a) 

' ' r p (H±i + 7 - a ) 



for ||C||<1, R e ( 7 -/})>^,Re(/?-a)>^. 

Substituting back we get (2.9), noting that one gamma is cancelled. 

3. RESIDUAL VARIABLES IN THE COMPLEX CASE 

Let X\ and X 2 be p x p hermitian positive definite matrix random variables defined in the complex 
field. Such a matrix Xj can be expressed in the following form: 

Xj = Xji + iXj2 

where i = V^T, Xji and Xj 2 are matrices with real elements such that Xji = Aj 2 > and Aj 2 = — Xj 2 . 
That is, Xji is symmetric positive definite and Xj 2 is skew symmetric. All the matrices appearing in this 
section are px p hermitian positive definite unless stated otherwise. 

If the growth and decay matrices in the complex field are respectively X\ and X 2 then the residual 
variable is 

Y = X 1 -X 2 . (3.1) 
Let Xi and X 2 be independently distributed matrix-variate complex gamma variables with the densities 

~9i(Xj) = |d ^ ( f j0 f J |det(^)l° f - p e- te ^^), X j = X;>0, Re( aj )>p-l (3.2) 

r p ( aj ) 

and gj(Xj) = elsewhere, where (•)* denotes the conjugate transpose of (•), Xj = X* > implies that Xj 
is hermitian positive definite, det(-) denotes the determinant of (•), |det(-)| gives the absolute value of the 
determinant and L p (-) is the matridx-variate gamma in the complex case, given by 

f p (a) = / |dct(A > )| Q - p e- tr ^)dA > (3.3) 

Jx=x*>o 

= 7r p ( p - 1 )/ 2 L(a)L(a - l)..T(a - p + 1), Re(a) > p - 1. 

Quadratic forms in complex Gaussian random vectors, giving rise to Whittakcr variables, are used in a wide 
variety of problems in communication and engineering areas, see for example, Hirasawa (1988), Divsalar et 
al.(1990), and Biyari and Lindsey (1991). 

Our aim here is to look at the residual variable in (3.1) when the input and output variables are 
independently distributed with the densities in (3.2). Then proceeding as in (1.7) to (1.10) we note that in 
order to compute the density of the residual variable one needs an integral to be evaluated, corresponding 
to the one in (1.10). This requires the definition of a Whittaker function of matrix argument in the complex 
case. Whittaker function of matrix argument in the complex case, denoted by W. t .(-), will be defined as that 
symmetric function, symmetric in the sense 

W a ,p(AB) = W a ,p(BA) 
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for p x p hermitian positive definite matrices, having the following integral representaton. 

f \det{Z)\ - a ~% |det(7 + Z)f +a S e~ tx{A2) dZ 

= \dct(A)\-P-i f p (/? - a + |) ei^W a , (A) (3.4) 

for Re(/3-a) > § - 1. 

As a consequence of the definition itself we can have the following result which will be stated as a 
theorem. 

Theorem 3.1. For A = A* > 

|det(/ + Z)r Q e- tr ^dZ 

Z=Z*>0 

= \det(Ap-e f p (p)ei^W^ M _ a+p) (A). 

One can establish a number of results on Whittaker function in the complex case which will be useful 
in working out the distributions of sums, differences and linear functions of matrix gamma variables in the 
complex case, evaluating probabilities associated with these distributions, evaluating the incomplete gamma, 
incomplete type-2 beta and other related integrals, and related problems. In order to illustrate the techniques 
one more result will be given here. 

Theorem 3.2. For B = B* > 0, U = U* > 0, M = M* > 0, Re(q) > f , 

|det(X + B)\ 2a - p \dct{X - U)\ 2q - p e- tT ^^dX 



L 



L 



x>u 

= f p (2q)|det(L> + B)\ a+q - p \dct(M)\- {a+q) 
xe^'l^W%_ 9)i(a+g _ f) (f), 
f = (U + B)iM(U + B)i. 



Proof. Consider the following transformations. 



Y = X-U, Z = (B + U)-*Y{B + C/)"5 
^dZ= \dct(B + U)\- p dY. 



Now the integral on the left side, denoted by /, is given by the following: 

f = e' t ^^\det(B + U)\ 2a+2q - p 

\dct(Z)\ 2q - p \dct(I + Z)\ 2a - p 



x 



Z=Z*>0 



xexp|-tr M(B + U)iz(B + U)i | dZ 



for Re(2q) > p — 1. Now interpreting with the help of the definition of a Whittaker function the result 
follows. 
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